import numpy
def miu(j,Y,pi,p,q):
    """
    j=0:n-1
    :param j:
    :param Y:
    :param pi:
    :param p:
    :param q:
    :return:
    """
    miu=(pi*(p**Y[j])*(1-p)**(1-Y[j]))/(pi*(p**Y[j])*(1-p)**(1-Y[j])+(1-pi)*(q**Y[j])*(1-q)**(1-Y[j]))
    return miu
def EM(Y,pi0,p0,q0,t):
    pi=pi0
    p=p0
    q=q0
    if t>666:
        return pi,p,q
    else:
        summiu=0
        summiuy=0
        sum1_miu=0
        sum1_miuy=0
        n=len(Y)
        for j in range(n):
            summiu=summiu+miu(j,Y,pi,p,q)
            summiuy=summiuy+miu(j,Y,pi,p,q)*Y[j]
            sum1_miu=sum1_miu+(1-miu(j,Y,pi,p,q))
            sum1_miuy=sum1_miuy+(1-miu(j,Y,pi,p,q))*Y[j]
        pi=sum1_miu/n
        p=summiuy/summiu
        q=sum1_miuy/sum1_miu
        t=t+1
        return EM(Y,pi,p,q,t)
Y=[1,1,0,1,0,0,1,0,1,1]
t=0
print(EM(Y,0.5,0.5,0.7,0))
